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ABSTRACT 

The self-consistent interaction between energetic particles and self-generated hydro- 
magnetic waves in a cosmic-ray pressure dominated plasma is considered. Using a 
three-dimensional hybrid MHD-kinetic code, which utilises a spherical harmonic ex- 
pansion of the Vlasov-Fokker-Planck equation, high resolution simulations of the mag- 
netic field growth including feedback on the cosmic rays are carried out. It is found 
that for shocks with high cosmic-ray acceleration efficiency, the magnetic fields be- 
come highly disorganised, resulting in near isotropic diffusion, independent of the 
initial orientation of the ambient magnetic field. The possibility of sub-Bohm diffu- 
sion is demonstrated for parallel shocks, while the diffusion coefficient approaches the 
Bohm limit from below for oblique shocks. This universal behaviour suggests that 
Bohm diffusion in the root mean squared field inferred from observation may provide 
a realistic estimate for the maximum energy acceleration timcscalc in young super- 
nova remnants. Although disordered, the magnetic field is not self-similar suggesting 
non-uniform energy dependent behaviour of the energetic particle transport in the 
precursor. Possible indirect radiative signatures of cosmic-ray driven magnetic field 
amplification arc discussed. 

Key viTords: acceleration of particles - (ISM:) cosmic rays - plasmas - instabilities. 



1 INTRODUCTION 

There is a growing wealth of observational evidence 
that non-linear amplification of magnetic fields occurs at 
the outer s hocks of supernova remnants, both upstream 



and down ( Achtcrbcrg ct ajj 19941: 
iBamba et al. 



Vink 
2011). 



L amind l2003l : 
This discovery 



i2005: Rako wski et al.l 
represents a significant development in the theories of parti- 
cle acceleration and the resulting non-thermal emission from 
supernova remnants. However, despite the fact that almost 
every attempt to model the non-thermal spectra from young 
supernova remnants in recent years has incorporated or in- 
ferred some level of field amplification, a detailed under- 
standing of the underlying mechanisms that produce such 
strong fields is still lacking. This can probably be attributed 
to the complexity of the problem, which, being highly non- 
linear, requires numerical modelling of processes occurring 
on a variety of length-scales. 

The amplification of magnetic fields to levels far in ex- 
cess of typical interstellar medium values is advantageous 
for a number of reasons. First, it seems essential to account 
for the x-ray synchrotron emission, its short time scale vari- 
ability, and its localisation to narrow filaments close to the 
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shock. Second, it is central to explaining the origin of cosmic 
rays with energies in excess of lO^'^eV. While the former can 
in principle be accounted for via fluid instabilities occurring 
at th e shock itself an d the region immediately downstream 
(e.g. iGuo et al.ll20il ). the latter requires strong fields up- 
stream of the shock that can only be produced from in- 
stabilities driven by energetic particles streaming ahead of 
the shock (e.g. Achtcrbcrg 1983; Boll 2004). Wh ile the non- 
reson ant hybrid (NRH) instabilities described in lBelll (|2004l . 
l2005h are believed to have the fastest growth rates, whether 
the resulting short-wavelength structures can significantly 
influence the highest energy particles has yet to be verified 
( BcU 2004; Rcvillc ct al. 2008). 

In this paper, we report on numerical simulations of the 
interaction between streaming cosmic rays and the thermal 
plasma upstream of the shock, including a self-consistent 
treatment of the cosmic rays. The thermal plasma is treated 
using an MHD fluid approach while the cosmic ray evolution 
is calculated using a spherical harmonic expansion of the 
Vlasov-Fokker-Planck equation. 

The outline of the paper is as follows. In the next sec- 
tion we introduce the relevant equations and the spherical 
harmonic expansion that is central to the approach used 
in this paper. From the resulting equations, the standard 
expressions for particle drifts and fluxes in the diffusion ap- 
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proximation are derived. These form the basis for many of 
the approximations required to self-consistently model the 
essential processes occurring in the precursor of an efficiently 
accelerating shock. Section 3 describes the numerical scheme 
and the simulation results are presented in section 4. We 
conclude with a summary of the main conclusions and ad- 
ditional discussion. 



in the local frame, in the sense that they are proportional to 
the angular component of the Laplacian in spherical coor- 
dinates. Since the spherical harmonics are themselves solu- 
tions to Laplace's equation, the resulting collision term takes 
a very simple form. 

Defining ff"^ = (/™)* the expansion takes the follow- 
ing form: Q 



2 PARTICLE TRANSPORT AND FLUXES 

The study of particle acceleration and magnetic field ampli- 
fication in supernova precursors is essentially one of magne- 
tohydrodynamics. Due to the large length-scales associated 
with the energetic cosmic rays, their interactions with the 
background field are mediated by the magnetic fluctuations 
supported by the background plasma. However, it is the cur- 
rents provided by the cosmic-rays that often play the domi- 
nant role in determining the evolution of the magnetic held 
in these precursors. We seek an effective numerical method 
to model the evolution of the cosmic rays in self-generated 
hydromagnetic waves and structures. For non-relativistic 
shocks Wsh <C c, it is generally accepted that efficient scat- 
tering maintains a quasi-isotropic phase spa ce distribution, 
at lea st to order Ush/c. In a previous paper (|Reville fc Belli 
|2012[ ). a particle in cell type numerical scheme was adopted 
to treat the cosmic ray evolution (see also [Zachar v 1987.; 
iLucek fc Belli I2OOOI ). Maintaining an isotropic distribution 
in multi-dimensions turned out to be computationally ex- 
pensive, with the number of particles per cell required to 
achieve the desired level of isotropy becoming prohibitively 
large for sim ulations in t hree s patial dimensions. 

Recentlv lBell et al.1 (|201ll ) performed test particle sim- 
ulations of cosmic-ray acceleration at oblique shocks in one 
dimension. The numerical scheme involved a spherical har- 
monic expansion of the distribution in momentum space 
based on the KALOS code previousl y developed for sim- 
ulations of laser plasma interactions l|Bell et al.l I2OO6I ). In 
this paper, we describe the full three dimensional expan- 
sion, which has been coupled to an MHD code. We adopt 
the mixed-frame coordinate system commonly used in as- 
trophysics, where momentum variables are measured in the 
local fluid frame 



[(p ■ V)u] . |i + ev ■ 
op 



B X 



dl 
dp 



= C{f) 



(1) 



accurate to second order in u/c l|Kirklll994l '). In this non- 
ine rtial frame an effective electric field term is introduced 
(cf. lSkillinelll975l ): 



E : 



du 
'dt 



+ {u ■ V)u 



We have also kept the third order term {u ■ V)u simply 
because the resulting expression for the electric field can 
be inferred from the momentum equation of magnetohydro- 
dynamics, which turns out to be convenient when included 
in the code. C{f) represent s the usual Fokker-Pl anck small 
angle scattering term (e.g. IChandrasekhaJ 1 19431 ) . The ad- 
vantage of using the mixed coordinates system is that the 
resulting collisions, as given by C(/), are elastic and isotropic 



f{p,e,^)^f2 J2 /r(p)pi™'(cos%'' 



(2) 



which, on substitution into equation Q, can be used to de- 
rive a set of coupled differential equations. These equations 
are quite lengthy, and for clarity they are written in single 
column form in the appendix. 

To motivate the various approximations that are made 
later in the simulations, we demonstrate how the well known 
transport equation is reproduced from the full expansion. 
For this purpose, it is necessary to retain only terms £ < 2. 
In the full numerical simulations, many more terms are used. 

Defining f = fo + fi ■ p/p where 

fo = fS, U = fi, fy = -2ml, h = -2Qfi 

and neglecting the terms associated with E above (which 
are second order in u/c), a somewhat tedious calculation 
shows that: 



m 



5 dxi ^ dp \p 



where ft = eB/'ymc is the directional rotational frequency, 
1/ is the collision frequency and summation over repeated 
index i — x,y,z is implied. In what follows we assume all 
particles are ultra-relativistic, and set everywhere v = c. 
Making the usual assumption that the first order terms react 
on a shorter timescale than the zeroth order component, we 
take a steady state solution for /i. Retaining only terms to 
first order in u/c leaves 

cV/o + n X /i = -ufi , 

which can be used to remove the /i dependence in the equa- 
tion ([3} to give 



dfo , l._ ^ dfo 



(5) 



+ V. 



Db 



[V/o + h{h- V/o) -hx V/o] 



where h — Cl/u is the hall parameter, and Db = c^/3v the 
usual Bohm diffusion coefficient. It can clearly be seen that 
in the relevant limits, the transport equation reduces to the 
well-known forms for parallel shocks 

2Il + u-Vfo = ^iV-u)p^ + V-iDBVfo) (6) 



^ For a mor e rigo rous discussion of this expansion see 
iTzoufras et al.l 1I2OI and references therein. 
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and perpendicular shocks 

^ + w ■ V/o = i(V ■ u)p^ + V-[D^{V-hx V)/o] 



where D± = Db{1 + h^y^ (|Bellll2008l ). 

We seek steady-state solutions, which will later be used 
to initialise the simulations. In the shock frame, the up- 
stream plasma moves in the positive x direction with uni- 
form speed Ush , and only gradients in this direction are con- 
sidered. Without loss of generality, a steady magnetic field 
is taken to lie in the x — z plane (i.e. h — hxX + h^z). The 
resulting upstream steady state solution is 



fo{x,p) = /o(0,p) exp / Kdx , with k 



(7) 



and the associated fluxes 



c 



-x + 



1 + hl' 



Mx,p) . (8) 



In the limit that the amplitude of the scattering waves 
(or u) does not increase significantly over a dista nce ^ , 
the re sults reduce to equations (5) and (6) of iBell et al] 
(|201ll ). However, small wavelengt h mo d es are in principle 
excited by the NRH instability of iBelll l|2004 ). with maxi- 
mum growth rate (for quasi-parallel shocks) 

1 Tier ^sh ^ 
7nr = n 

2 nth VA 

such that the number of exponential growth factors in am- 
plitude of magnetic fluctuations over a distance equal to the 
fluid crossing time of is 



hl + hi 



nth Ush VA 6 1 + h'^ 

For efficient shock acceleration, this number can exceed 
unity. In fact, at the highest energies, it must exceed unity 
if magnetic fields are to be amplified upstream. 



3 NUMERICAL SCHEME 

The numerical scheme is a natural extension of earlier hy- 
brid MHD-k i netic w ork based on the MHD code used in 
Reville et al.l 120081), with the particle-in-cell approach of 
ReviUe fc Belil l|2012l ) replaced by a Vlasov-Fokker-Planck 
routine. The full expansion of equation (lAip is solved nu- 
merically using a third order Runge-Kutta method. Making 
use of the orthogonality relations of the spherical harmonics, 
it follows that the cosmic-ray particle and current density, 
measured in the local fluid frame, are 



n^q d^pf = 47V / dp//o° 



and 



d^pvf ■- 



—q j dpp V 



f? 
29(/i') 



(9) 



(10) 



The MHD equations, modified by the pre sence of cosmic 
rays have been described in detail elsewhere (|Zacharvll 19871 : 
[Beiill2004l ). but we repeat them here for completeness. The 
standard MHD momentum equation is supplemented with 



the additional effect of the return current induced by the 
cosmic-ray streaming, giving 

du 



dt 



4tv 



B X (V X B) j„ X B 



Defining the MHD energy density e = pu^ /2 + P/(7 — 1) + 
B'^/S-K, where 7 is the adiabatic index of the background 
gas, it follows that the energy density also differs from its 
usual conservative form 



de 
di 



+ V. 



e+p + 



Stt 



{u B)B 



E (11) 



i.e. if the background plasma can set up an electric field to 
oppose the cosmic ray current, work is done by the cosmic 
rays on the MHD fluid. Note that there is no need to trans- 
form back to the laboratory frame since the return current 
is always defined locally to oppose that of the cosmic rays. 

Since the primary focus of this paper is the growth of 
magnetic waves and its feedback on the particles driving the 
growth, we make the following approximations: 

(i) To reduce memory requirements and allow for a larger 
spatial domain as well as larger number of harmonics, we 
dispense with the momentum g rid. This is achieved apply- 
ing the same technique used in iBell et al] (|201ll ) , where a 
uniform power-law distribution is assumed for all harmonics 



9/; 



rm 

1— 

P 



(12) 



This is similar to making a mono-energetic approximation, 
but ensures the Ex B drifts remain properly accounted for. 
For the parameters we consider, this is the most important 
effect. To prevent divergence of the cosmic-ray pressure, we 
have taken the value s = 4.1 in all simulations. This also 
ensures all ff™ terms (see appendix, equation IA5|) do not 
vanish exactly. 

(ii) a finite collision frequency !/(<g SI) is chosen in order 
to guarantee the closure of the spherical harmonic expan- 
sion. As can be seen from the appendix (Eq. IA6|) . in the 
absence of other effects the high harmonics are damped ex- 
ponentially 

fTit + At) = frit) exp [-'^£{1 + l)At 

It is therefore possible to truncate the harmonic expansion 
at ^ = imax, subject to the condition that Lmax > ^^fl/u. It 
will be demonstrated that after the field has been allowed to 
grow sufficiently, the self-consistent scattering on magnetic 
field structures dominates the artificial collisions imposed by 

V. 

(iii) Perhaps the most crucial approximation that we 
adopt, which is a theoretical tool quite unique to the spher- 
ical harmonic expansion technique, is to add an additional 
term proportional to the large scale gradient of the isotropic 
part of the spectrum fhsix). This results in an additional 
source term in the ^ expression for fi f Equation ! A2I) . i.e. 
we include 



A" 



dx 



This term acts like an external driving force, which allows us 
to model the larger scale dynamics of the precursor. In the 
absence of self- generated magnetic fields, this term would 
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allow the system to relax to the standard steady state solu- 
tion, where advection and diffusion balance. In order to use 
this approximation, the gradient must vary on scales larger 
than the length of the simulation domain Laim, ie. 




Similarly, the use periodic boundary conditions for all other 
quantities demands that this condition be satisfied, since 
otherwise the large scale gradient is resolved, and different 
boundary conditions or full shock simulations are necessary. 
The gradient that we use in the simulations is taken to be 
that associated with the steady state solution described in 
the previous section, where the value of Db is determined 
from the numerical value for v used in the simulation. 

Previous PIC/Hybrid kinetic simulations of cosmic-ray 
streaming instabilities including feedback on the driving par- 
ticles, have suffered from the fact that the energy must be 
transferred to the back ground plasma at th e expense of the 
cosrn i c-ray anisotrop y jLucek fc Belli I2OO0I : IStroman etaP 
I2OO9I : iReville fc Belli I2OI2I '). Hence"~onlv a finite amount of 
free energy, determined by the initial conditions , was avail - 
able in the system to amplify the mag netic field. Hii (j2004l l 
circumvented this problem by noting that on scales much 
less than the gyroradius of the cosmic rays, the feedback on 
the particles was, at least in the linear theory, negligible, 
and a constant uniform current could be used. In this paper 
we bridge the gap between these two regimes, by using the 
local approximation above, which allows us to control the 
driving of the cosmic ray current. 



4 SIMULATION RESULTS 

Despite the fact that a spherical shock expanding into a 
large scale uniform magnetic field has considerably more 
of its shock surface at large obliquities, parallel or quasi- 
parallel shocks have received by far the most attention. 
This may be largely due to the belief that injection into 
the acceleration process is suppre ssed at high Mach number 
quasi-perpendicular shocks (e.g. iBaring et al.|[l993 ). How- 
ever, observations of supernovae would suggest that this 
is not the case. Shell-type supernovae are quite abundant 
in our galaxy, and it is extremely unlikely that these rims 
represent exclusively quasi-parallel shocks. Kinetic simu- 
lations of th e shock m icrophysics are provi ding new in- 
sight (see e.g. lGargate fc Spitkovskv 2012; Mats umoto et all 
|2012| . and references therein), but the solution to the injec- 
tion problem in its entirety remains unsolved. In the follow- 
ing, we assume that the shock has already put a considerable 
amount its energy into the cosmic-ray population, but not 
yet significantly amplified the magnetic fiuctuations. 

The simulations are performed in three dimensions, us- 
ing periodic boundary conditions, where the computational 
domain represents the evolution of a small region at rest in 
the precursor of a supernova remnant shock. Using the lo- 
cal approximation described in the previous section, we in- 
vestigate the effect of the angle between the shock normal, 
as represented by the large scale gradient, and the mean 
magnetic field. A turbulent component, {SB^)/Bq « 1%, is 
added to the background magnetic field, to seed the growth 



of waves. This seeding is essential as VFP codes do not suf- 
fer from noise. All simulations begin with the plasma com- 
pletely at rest, and plasma beta of unity. While the simula- 
tions attempt to model as close as possible the parameters 
relevant for young supernovae, the separation of scales is al- 
ways a challenge. This manifests itself through the different 
timescales for the MHD updates, determined by the maxi- 
mum fast-magnetosonic velocity, and the cosmic-ray veloc- 
ity, which is essentially the speed of light. Thus the CFL 
condition on the VFP equation requires a certain amount 
of sub-cycling. To allow the simulations to run in a rea- 
sonable time we are forced to use an initial Alfven velocity 
somewhat larger than expected in a real situation, typically 
< lO^'^c. Thus, most of the simulations consider shocks with 
Alfven Mach numbers on the order of Ma ~ 50 — 100. The 
code is terminated when the maximum fast magnetosonic 
speed approaches the speed of light, which usually results 
due to rapid non-adiabatic heating in the non-linear regime. 
This non-adiabatic heating may have observational signa- 
tures which will be discussed in section [431 

4.1 Oblique shocks 

The velocity of the intersection point between the shock and 
a given magnetic field line is Uint = Msh/ cos 6m, where 9n 
is the angle between the large scale field and the shock nor- 
mal. If the fields are strongly amplified to values SB /Bo > 1, 
this quantity is, at least locally, quite poorly defined. How- 
ever, since we start with quite modest fiuctuations, we refer 
to the obliquity with regards initial mean field. Already at 
Ush = c/10 an angle 9m > 85° is required for a superluminal 
intersection velocity, such that even at quite large obliqui- 
ties, particles should be able to escape upstream along field- 
lines. 

Starting with the smaller obliquities, figure[l](a) shows 
the time evolution of the fluxes and turbulent fleld compo- 
nent for a shock at an initial angle 6m ~ 60° to the shock 
normal. The cosmic-ray fluxes grow quickly from zero to 
match the values predicted from Eq. ([S]). In fact the peak 
values correspond almost exactly to these values, suggesting 
that the diffusion approximation is reasonable at this point, 
and that the scattering is not yet dominated by the self 
generated fields. This is to be expected since fiux-freezing 
determines that the fields cannot grow before the plasma 
is set in motion, so that th ere is an i nitial inertial phas e 
where no growth is observed (|Bellll2004l : iReville et allbOOSl I , 
and the cosmic-ray scattering is likewise unaffected in the 
weak fiuctuations. However, once the field starts to grow, 
the magnitudes of all components of the cosmic-ray current 
begin to decrease. By the time the field becomes completely 
disorganised SB ~ Bo, the shock begins to behave gradually 
more like a parallel shock, for which the parallel component 
of the current dominates. This parallel component continues 
to grow on account of the strong gradient which assumes a 
fixed value for k (equation [7| relevant for a 60° angle be- 
tween the shock normal and the magnetic field. 

Increasing the magnetic field angle to 80°, it can be 
seen that the time evolution of the various components 
of the current changes slightly. Figure [1] (b). The parallel 
component is found to grow slowly at first, but then over- 
shoots the value predicted by the diffusion approximation, 
Eq. ((Sj. The overshoot corresponds approximately to the 
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fluctuation level where the self-consistent scattering domi- 
nates (SB /Bo > v/Q), at which point the cosmic rays dif- 
fuse across the field li nes more easily . Prev ious VFP sim- 
ulations performed bv iHornsbv et all (|2010l ). investigating 
cross field transport in the presence of synthetic turbulence, 
found that cross field transport approaches the Bohm limit 
at SB ~ 2Bo. Returning to the results shown in Figure [1] 
(b), it is found that the value for jy also overshoots, while 
jz falls short of its predicted value. This is also consistent 
with a very rapid increase in the collision frequency. This 
may be due in part to the zero-th order ja x Bo force ac- 
celerating the plasma almost instantly as the angle grows. 
The currents appear to saturate in a very similar fashion to 
the previous run at 60° (compare Figures [1] (a) and[T](b)). 

A similar behaviour occurs at the extreme limit of su- 
perluminal shock propagation, 9n = 90°, where both the x 
and y components of the current overshoot, with jx gradu- 
ally increasing into the non-linear regime before saturating, 
while jy decreases towards zero, Figure [1] (c). Again, the 
driving term is fixed in time so the gradient is very steep. 
The regulation of this gradient is discussed further in section 

In all cases, the final stage of evolution appears to be 
similar, namely that cosmic rays are prevented from slid- 
ing along the field lines, as represented by jz, while the 
drifts orthogonal to the plane of the mean field are gradually 
damped. This suggests a universal behaviour for efficiently 
accelerating shocks, such that they tend asymptotically to- 
wards parallel-like shocks. While this phenomenon has been 
widely believed to occur, this is the first time that it has 
been demonstrated in a self-consistent manner. The global 
effect this has on the acceleration is uncertain at present. 

A possible by-product of this behaviour may reveal itself 
on interaction with the shock surface itself. Previous numer- 
ical simulations by Zirakashvili & Ptuskin (2008) demon- 
strated that for a parallel shock, the non-linear structures 
generated by the NRH instability can, on crossing the shock, 
result in stretching of the magnetic component parallel to 
the shock normal, such that it dominates over the compo- 
nent in the shock plane. We find that similar structures are 
found in simulations at all obliquities, suggesting that even 
initially perpendicular shocks may result in radially point- 
ing magnetic fields downstream of young supernov a rem- 
nant shock, as has been observed in several remnants (|Milnel 
ll98'iD . 

A caveat regarding the non-linear evolution concerns 
the finite lifetime of the simulation region, being in the 
upstream rest-frame, it must at some stage be advected 
through the shock. This problem is most severe in the case 
of highly oblique shocks, where the precursor is initially on 
the order of only a few gyroradii. However, as the magnetic 
field becomes increasingly disordered, particles will pene- 
trate deeper into the upstream plasma, and allow further 
time for magnetic field growth. Ultimately, for shocks ef- 
ficiently accelerating cosmic-rays, self-generated fields take 
over and the shock will behave more like a parallel shock. 
However, sufficiently far upstream, the fields will remain reg- 
ular, at least at the level of pre-existing interstellar turbu- 
lence, and particle escape to infinity is impeded. Clearly full 
shock simulations are required to treat this problem self- 
consistently. This will be addressed in a follow-up paper. 

It is tempting to connect these results to the previous 




(c) 6»jv = 90° 










Jz 















' time 

Figure 1. The lab frame current densities j = ricrsu + j' and 
growth of magnetic field as a function of time. Initial conditions 
are chosen such that u = O.lfi, iigh = 0.03c. Magnetic field is in 
units of Bo. Time units are 20.8/Qo i-e. if the CR are at 100 TeV, 
and So = 3 /iG, time is approximately in years. 



work of iBell et al.l (|201ll ). where it was demonstrated that 
the shape of the non-thermal power law could deviate from 
the standard test-particle value at oblique shocks, depend- 
ing on the angle and collision frequency. This deviation is 
essentially a direct consequence of the inability to match 
the perpendicular currents across the shock in the diffusion 
approximation, which demands fo{x > 0,p) = const. We 
have just demonstrated that the time-asymptotic state of 
any efficiently accelerating shock front is one in which the 
perpendicular currents vanish, which implies the existence 
of a possible radiative signature that can be used to probe 
the level of upstream field amplification. To illustrate this 
point, SN 1987a is taken as an example. 

Considering very young supernovae, the results pre- 
sented here suggest that shocks which are efficiently accel- 





Figure 2. Slice through the x — y plane demonstrating the correlation between Ax and ricr- While this is found to occur at early times, 
no additional focussing was found in any of the simulations, and the correlation eventually disappears in the non-linear stage. 



sive suite of simulations investigating the growth of magnetic 
field in this configuration, to investigate a number of effects. 

Recently, iReville fc Belil (120121 ') presented a theoretical 
and numerical investigation of cosmic-ray filamentation and 
its role in the generation of large scale magnetic fluctuations. 
The results were based on a two-dimensional analysis that 
assumed slab-symmetry along the direction of cosmic-ray 
anisotropy. It was demonstrated that the cosmic-ray den- 
sity correlated with the component of the magnetic vector 
potential parallel to the direction of cosmic ray streaming. 

Using the VFP code, we are able to drop the approxima- 
tion of slab symmetry and investigate if the process occurs 
when the third dimension is included. The code is set up as 
before, with periodic boundary conditions and a large scale 
imposed gradient to drive a current consistent with equa- 
tion ([8}. It is found in all cases that the cosmic- ray density 
correlates with the vector potential during the linear growth 
{5B <^ Bo) as can be seen in Figure (2] In fact, the correla- 
tion also exists at oblique shocks, although the correlation 
is between the cosmic-ray density and the component of the 
vector potential parallel to the mean cosmic-ray curren t. In 
this sense, the linear analysis of iReville fc Bel] l|2012t ) still 
applies, and one might expect growth of large scale fields. 
However, the resulting filaments that form do not continue 
to focus, and the structures are ultimately destroyed in the 
non-linear stage. It should be pointed out that we only con- 
sider the evolution of cosmic rays at a single energy, and that 
the simulation domain is periodic. It is still possible, and in 
fact quite likely, that the leading ed ge of the precursor can 
still fi lament (see recent results from lCaprioli fc Spitkovskvl 
(120131 )'). The periodic simulations presented here are more 
representative of a region deep in the precursor. An exten- 
sive parameter scan failed to produce noticeable elongated 



erating cosmic rays might have an observable time depen- 
dent behaviour, as both the cosmic-ray efficiency and self- 
generated magnetic field fluctuations increase. SN1987a rep- 
resents an interesting example, since the radio flux has been 
increasing steadily for over 20 years, implying an increasing 
acceleratio n efflciency, while th e spectrum has been steadily 
flattening (|Zanardo et al.|[2010l ). Reconciling this behaviour 
within the non-linear diffusive shock acceleration framework 
requires a somewhat arbitrary flne-tuning of the electron- 
proton injection ratio (|Berezhko et al.ll201ll ). However, this 
time-dependent behaviour is quite consistent with accelera- 
tion at an oblique shock. Since t he shock is believed to b e 
expanding into a Parker spiral (|Kirk fc WassmannI Il992l ) , 
at least for some of its lifetime, the majority of the shock 
surface should have large magnetic obliquities. As the accel- 
eration becomes more efficient, the magnetic fluctuations be- 
come amplifled, and t he collision f requen cy increases. Com- 
paring to flgure 1 of iBell et al] (|201ll ). this will lead to 
a gradual flattening of the spectrum. As the held enters 
the non-linear regime of magnetic field amplification, the 
perpendicular currents are damped near the shock, and a 
power-law spectrum closer to the standard test-particle re- 
sult should result. 'While this explanation is clearly not in 
any way rigorous, it does provide an alternative to the non- 
linear shock acceleration model, which neglects the effects 
of the obliquity of the upstream magnetic fields, and addi- 
tionally any deviations from the diffusion approximation. 

4.2 Parallel shocks 

Parallel shocks have received the greatest attention in the lit- 
erature, both in terms of particle acceleration and resulting 
magnetic field amplification. We have performed an exten- 
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structures in our simulations. On the other hand, the re- 
sults do seem to suggest that cosmic rays can only escape 
the source at the highest energies, which is consistent with 
current models of cosmic-ray escape. Full, energy- dependent 
shock simulations are still required to rule out the possi- 
ble significance of cosmic-ray filamentation as a low energy 
cosmic-ray escape channel. 

The simulations allow us also to address, the interaction 
between the cosmic rays and the self-consistently amplified 
magnetic fields, over a spatial range that covers both the 
gyro-resonant interactions and the sub-Larmor scale fluctu- 
ations where the fast growing NRH instability operates. A 
number of interesting effects are found. First, it is noted that 
the evolution of the plasma is quite different from that typ- 
ically seen in the fixed c o smic-ray current simula tions (jBelll 
120041 : iReviUe et ai]l2008l : IZirakashviU et al.ll2008l ). Figure [3] 
shows the evolution of the different components of the en- 
ergy. Previous MHD simulations of the field amplification 
using a constant uniform current typically go through three 
stages. An initial inertial phase, where the plasma is set in 
motion, and the magnetic field is either constant or can ac- 
tually decrease as energy is transferred to the background 
fiuid motions. Once in motion, the plasma starts to stretch 
the field lines, and reinforced by the j x B force, leads to 
a run-away instability. Finally, in the non-linear phase, the 
plasma motion is highly disordered and rota tional. The ki- 
netic energy continues to grow as fiPelletier et al.ll2006l : 
IZirakashvili et al.ll2008l ). and the magnetic energy grows in 
equipartition with the same scaling. 

When we include the feedback on cosmic rays, the initial 
behaviour is the same, but the growth of kinetic energy stag- 
nates at early times. The magnetic field starts to grow, with 
total energy density in the growing magnetic modes quickly 
dominating the kinetic energy. Already at SB^ /Bq ~ 0.1 the 
scattering on magnetic structures acts to reduce the current 
such that the driving term can no longer compensate. Even- 
tually, the cosmic-ray current becomes too weak to have any 
infiuence on the magnetic field on a timescale comparable 
with the precursor crossing time. 

An interesting result of amplifying the magnetic field 
on sub-Larmor scales, is that the frequent scattering on the 
resulting non-linear fiuctuations prevents the current from 
following the field on large scales. This can clearly be seen 
in Figure |4l where the dominant magnetic field structures 
are on scales much less than the particle gyroradius which 
is 3/4 the box length and 3 times the box width. Thus, one 
can safely assume that the usual expression for the linear 
growth rate of parallel modes (fe||i?o) 



7 = 



kjciBp 
pc 



(13) 



holds up to wavelengths larger than the gyroradius of the 
particles driving the growth. Of course, the simulations are 
initialised with a perfectly uniform mean field, and varia- 
tions in the orientation of the large scale magnetic field 
in the pre-existing ambient circumstellar medium may in- 
troduce interesting additional effects (Giacinti ct al. 2012). 
Nevertheless, a faster growth rate at long wavelengths is 
clearly implied. 

The role of c ollisio nal effects has also been discussed 
m ISchure fc Belil i 201 ll ). Whether this is the same insta- 




time 



Figure 3. Time evolution of the different components of the en- 
ergy density. Usb = — ^o)/^''^ energy density in the 
magnetic fluctuations having subtracted off the energy associated 
with the mean field. All energy densities are in units Bq/Stt. (The 
sharp jump in jx from t=0 to t=l is due to the fact that the cur- 
rent is only stored at integer timesteps.) 



bility is uncertain, since non-resonant scattering on short 
wavelength fluctua tions appears to dom inate, and unlike the 
modes analysed in lSchure fc Belli ()201l| ). the growing waves 
are exclusively of a single polarisation, rotating in the op- 
posite sense to the cosmic rays in the mean fleld, i.e. non- 
resonant. 



4.3 Regulated driving 

For the simulations reported in the previous sections, a fixed 
driving term, determined by the initial conditions and our 
numerical collision term was used (equation [71 . Since the 
purpose of these investigations is to study the interaction of 
particles with their own self-generated non-linear structures, 
we naturally expect a reduction in the mean free path, and 
consequently, for the system we consider, a similar reduction 
in the current (Figure [3]). 

At an actual (quasi-parallel) shock, the increased scat- 
tering rate would prevent particles escaping far upstream of 
the shock, and the scaleheight of the precursor would de- 
crease. Thus, it is reasonable to expect that the cosmic-ray 
gradient and the diffusion might conspire to maintain the 
total current density at a constant level, i.e. 



Jcr 



d pD — 

ox 



constant 



We have performed simulations where the large scale gra- 
dient (the driving term) is increased in such a way as to 
maintain an approximately constant cosmic-ray current as 
the field continues to grow. Starting from the usual result, 
we allow the current to build up to t = 3, and store this 
value as a reference current. As the magnetic fiuctuations 
grow and dominate, the driving term is increased using a 
normalisation parameter A. This parameter then becomes a 
measure of the effective reduction in the diffusion coefficient 
A — Db/D, where Db is the Bohm value determined from 
the initial conditions, which for the simulations shown has 
= $1/8. It is found that, not only does the diffusion co- 
efficient reduce well below the value imposed by the initial 
conditions (as expected), but by the end of the simulation. 
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Figure 4. Plot of the field lines for j and B from time t = 10 for the simulation shown in Figure [S] The approximately straight dark 
lines correspond to the current, while the green helical lines show the magnetic fieldlines. The gyroradius of the cosmic rays in the mean 
field is 3/4 the length of the simulation box. 



has even reduced below the Bohm limit in the pre-existing 
field VbS ~ 20. This can be clearly seen in Figure[6l The sim- 
ulation is terminated when the maximum fast-magnetosonic 
velocity becomes comparable with the cell crossing time of 
the cosmic rays. This happens quite early, due to the artifi- 
cially large Alfven velocity used in the simulations. In reality, 
the fields would conti nue to grow for considerably longer, as 
found for example in iBelj (|2004l 'l , and the diffusion coeffi- 
cient would likely continue to decrease in parallel. We note 
however, that the scattering is likely dominated by localised 
regions of very intense magnetic field. At the end of the sim- 
ulation the magnetic structures are sufficiently strong and 
on sufficient scale that they can deflect a cosmic ray through 
~ 90°, as shown in Figure (5] 

Since the maximum value of the magnetic field is in 
general time-limited (saturation of the magnetic field has 
not been found in any simulation with constant /sustained 
current), it is expected that the diffusion coefficient of the 
cosmic-rays will continue to decrease with growing magnetic 
field, presumably well beyond the level shown here. It is com- 
monly assumed that Bohm diffusion in the amplified mag- 
netic field is an acceptable value for the diffusion coefficient 
in most calculations of relevance to diffusive shock accelera- 
tion. Thus, it is important to be clear about what this means. 
In section [21 the Bohm scaling argument was derived. Typ- 
ically when referring to the Bohm limit, it is thought that 
the collision frequency matches the gyrofrequency, such that 



For highly disordered magnetic fields, confusion can now 
arise, as it is unclear what value of the magnetic field should 
be taken when determining the gyroradius, and more im- 
portantly, if the Bohm limit itself applies in nature in the 




20 40 60 80 100 120 



X-axis 

Figure 5. Slice through the x—y plane, showing the magnitude of 
the magnetic field, in units of Bq. The regions of intense magnetic 
field have an average strength of ~ 25i?o and a size L ~ 10. The 
gyroradius of the cosmic rays in the initial field is Tgo = 256, such 
that these regions of intense field can deflect particles through 
Ae ~ L/r^ ~ 90°. 



presence of strong amplification. In the simulations pre- 
sented here, our experimentally determined diffusion coef- 
ficient is half that of the Bohm limit value if we adopt the 
gyroradius in the pre-amplified magnetic field, but roughly 
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a factor of four times the Bohm limit in the root mean 
squared field0 However, the diffusion coefficient was still 
reducing at the time the simulation was terminated, sug- 
gesting that we might expect the diffusion coefficient at a 
supernova to correspond to its Bohm value in the amplified 
magnetic field, modulo a numerical factor of order unity. 
Thus, it may not be such an extreme approximation when 
considering the maximum energies at young supernovae, 
to adopt the Bohm limit in the amplified magnetic field, 
namely the one inferred from observations. Unfortunately, 
it remains uncertain as to whether the majority of the am- 
plification occurs in the upstream, or acts simply as a seed 
for Ric htmeyer-Meshkov type instabilities at the shock fron t 
itself (jGiacalone fc Jokipii 2007 : Schure fc Belli lin prep.l ). 
However, as noted by IGuo et al. (|2012l '). the growth in the 
downstream region due to such clumpy interaction occurs 
over a large distance, and some additional mechanism is 
required to explain the strong fields observed close to the 
shock. This favours a scenario in which at least a sizeable 
fraction of the amplification occurs upstream. Future sim- 
ulations that include a shock in the simulations will shed 
further light on this issue. 

Finally we note that young supernovae having clear 
Balmer emission from the outer shock may be used as an 
indirect probe of upstream amplification. All of the simu- 
lations, independent of magnetic obliquity show a common 
feature, namely that in the non-linear regime, the thermal 
energy density is greater than the magnetic energy density 
by typically an order of magnitude, e.g. Figure |6l Thus, if 
fields of 100/xG are produced upstream, a thermal energy 
density in excess of f keV/cc may also be produced. The sur- 
vival of neutrals in such high temperature plasmas is highly 
unlikely, unless the precursor is very short. However, a short 
precursor will not provide sufficient time to amplify fields 
to such high values. Several young remnants show regions 
of both non-thermal x-ray filaments, synonymous with effi- 
cient particle acceleration, and broad/narrow Ha emission, 
suggesting survival of neutrals into the downstream (e.g. 
Winkler fc Long||l997l : iGhavamian et al.ll200ll : iHelder et al] 



20091 ). If indeed the magnetic fields are produced predom- 



inantly upstream, these regions should anti-correlate, thus 
providing indirect evidence of in-situ acceleration of cosmic 
ray protons or nuclei, or indeed the lack thereof. Further 
observations are necessary to answer this question. 



5 DISCUSSION 

We have presented a series of numerical simulations inves- 
tigating the interaction between cosmic rays and magnetic 
fields in the precursors of supernova remnants. The inves- 
tigations are carried out using a hybrid MHD-kinetic code 
which applies a spherical harmonic expansion of the Vlasov- 
Fokker-Planck equation. This code allows us to consider a 
range of multi-scale phenomena previously inaccessible to 



^ This applies only to the particles driving the field growth 
within our mono-energetic approximation. Previous work by 




iReville et al.l | |2008| ) have already demonstrated the energy de- 
pendence of the transport properties in self generated turbulence 
when the particles were all at considerably lower energies than 
those assumed driving the field growth. 



time 



Figure 6. Growth of the different components of the energy den- 
sity as a function of time for simulations with regulated driving. 
The current density jx and A the renormalisation parameter of 
the driving are also shown, i/ = 0.125r2o, u^i^ = 0.1c. Time units 
are 5.2/Oo 



numerical simulations. A number of results have followed 
from these simulations that can be applied to existing the- 
ories of particle acceleration, and in particular the observa- 
tional signatures of efficient acceleration at young supernova 
remnants. 

Previous work by iReviUe fc Bel (|20f2l ) proposed a 
method whereby low energy cosmic rays could escape the 
remnant via self-focussing into low density cavities sur- 
rounded by strong magnetic fields. The a nalysis and simu- 
lations performed in iReville fc BeUl l|2012l ) were exclusively 
in a 2D slab geometry, which limited the amount of focus- 
ing that could be achieved. This is a consequence of the 
conservation of the component of the vector potential into 
the plane in such a geometry. Including the third dimen- 
sion allows for growth of this component, and in theory, 
continued focusing, which could result in formation of fila- 
mentary structures aligned along the direction of cosmic-ray 
streaming. The curr ent 3D s imulations appear to support 
the linear analysis of lRevill~fc Bell ( 201 2i) . but that the re- 
sulting structures become disrupted over a short distance. 
While this disfavours the possibility of non-diffusive escape 
channels for cosmic rays from the remnants, it does pro- 
vide convincing evidence that the escape of cosmic rays into 
the interstellar medium at quasi-parallel shocks is limited 
to a narrow window around the maximum energy, consis- 
tent with current models. 

Oblique and perpendicular shocks present a slightly dif- 
ferent view, although there is clear evidence for a universal 
behaviour closer to the shock. Far upstream of the shock, 
where the scattering is weak, particles are tied to field-lines, 
and advected back towards the shock. It is demonstrated 
however, that the self-generated fields closer to the shock 
act to damp any diamagnetic-currents and in a sense, behave 
very much like parallel shocks. At the same time, consider- 
able amplification of the magnetic field can occur. Previous 
suggestions that quasi-perpendicular shocks are faster ac- 
celerat ors have imp licitly assumed that the scattering was 
weak (|jokipiilll987l ). Later calculations demonstrated that 
the maxim um energy became a p roblem of geometry rather 
than time (|Achterberell2004l : iBelil 12008). as the energisation 
relies on drifting along the shock surface. On the global scale 
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of the precursor, we now have a rather unusual situation, 
quite distinct from parallel shocks, where the particles can 
diffuse isotropically close to the shock, but are strongly con- 
fined by the upstream undisturbed field. Whether this acts 
to increase or decrease the maximum energy is not imme- 
diately clear. However, full shock simulations will be per- 
formed to answer this question in the near future. 

Oblique shocks may also play an important role in ex- 
plaining observations from young supernovae, or likewise, 
early observations can provide important information on 
the acceleration process. Immediately after break-out, su- 
pernovae typically expand into a progenitor wind. It is usu- 
ally assumed in such winds that the radial magnetic field 
falls off quite rapidly, while the perpendicular component 
will persist. Thus young supernovae provide ideal labora- 
tories to investigate many of the ideas presented here. We 
have demonstrated, albeit in a very qualitative sense, how 
the current evolution of SN 1987a can be explained within 
the framework of acceleration at highly oblique shocks with- 
out the need for non-linear cosmic-ray pressure feedback. 

The possibility of sub-Bohm diffusion (with respect to 
the pre-existing ambient field strength) of cosmic rays at the 
highest energies has been demonstrated self-consistently for 
parallel shocks, under the assumption that the driving reg- 
ulates itself. On the other hand, the normal component of 
the diffusion tensor at oblique shocks is initially sub-Bohm, 
but tends to increase due to amplification of magnetic per- 
turbations. Regulation is more complicated in these geome- 
tries, and full shock simulations are required. However, the 
high resolution simulations presented in this paper, show 
that the self-generated fields naturally produce highly dis- 
ordered field lines strongly suggesting a universal behaviour 
of isotropic diffusion at shocks that are efficiently accelerat- 
ing cosmic rays. 
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APPENDIX A: SPHERICAL HARMONIC EXPANSION OF VLASOV-FOKKER-PLANCK EQUATION 

In the following, we choose the polar axis to lie along the s-direction, such that the Cartesian momentum vector expressed in 
terms of spherical angles is p = (pcos ^, psin S cos (;/>,psin ^ sin 0). Using the recursion relations and orthogonality condition of 
the spherical harmonics, the time derivative for each component of the expansion can be written in the form 



dfl 
dt 



b: 



E. 



(Al) 



Most of these terms have been previously given in iBell et al.1 (|2006l ) and lXzoufras et al.l l|201ll '). but are repeated here so 
all the equations can be found in one place. The mixed coordinate frame introduces several additional terms, contained in 
the T^p i, which can be associated with the adiabatic term in equation ^ in the main text. The subscripts a,l3 £ {x,y,z} 
correspond to the relevant term containing the partial derivative dua/dxp. Although cumbersome, the derivation of each 
term is straightforward, and the final expression often displays a noticeable symmetry, a clear indication of the applicability 
of spherical harmonics in solving problems of this nature. 

The advection terms are as follows 
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where D± = |- ±i|-. 

oy oz 

Defining the cyclotron rotation vector 17 = qB/'ymc, with q the particle's charge and 7m its (relativistic) mass, the 
rotational terms associated with the magnetic field components are 



sr^" = in,mfr + ^ [(a+if^y)fr"' 

B'-I = -l{£ + 1)3? [(fi, - iny)fl] 



(£ + m m)(f7z - if^y)ff 



(A3) 



The spherical harmonic expansion has the additional feature that th e rotational ter ms are algebraic, and do not involve partial 
derivatives, as would be the case for a Cartesian tensor expansion (|johnstonlll"960h . 

Recalling that in the mixed coordinate frame, the local acceleration has an equivalent form to the electric field in the 
standard Vlasov equation, we thus use the notation 
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from which it follows 
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where 



-U+l) " / {^+1) rm 



/r) 



(A5) 



Since the adiabatic term has two occurences of the momentum vector, it contains nine different terms associated with the 
different combinations of the axes. However, due to the symmetry about the polar axis, most of these t erms are very similar 
and can be grouped together. Some of these terms have been previously derived in one (|Bell et al.ll201ll ) and two dimensions 
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(C. Ridgers private communication). 
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dy dz J [ (2£+l)(2£ + 3) 



(2£- 1)(2£+ 1) 
■ m)(£ - m - 1)(£ - m - 2)(£ - m - 3) 



(2£-3)(2£-l) 



^m+2 
^f-2 



l)+m2 



(2£+ 1)(2£- 1) 



HV 



(£ - m)(£ - m - 1) (£ + m + 1)(£ + m + 2) ' 

-(jr^_2 fo/i ■ rr\/o/; ■ o\ -"^+2 



(2£- l)(2£-3) 



(2£ + 5)(2£ + 3) 



dy J 2£ + 



where 
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For m < 2 we have 
iii + 1) 



■ pY 



(2£- 1)(2£+ 1) 



(Hir 



* and 

{i + 2){i + 3) 
(2^ + 3)(2£ + 5) 



(2^+l)(2£ + 3) 



(2£-3)(2£-l) 



{Gur 



Finally, since the spherical harmonics are themselves solutions to Laplace's equation, the collision term takes a particularly 

simple form 



(A6) 
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